ATM 623: Climate Modeling

Lecture 2: Solving the zero-dimensional energy balance model

A self-study exercise in lieu of lecture

About this document:

This is an IPython Notebook. It combines formatted text, Python code, and output generated by the code (including graphics).

For an introduction to the Notebook format, see http://ipython.org/ipython-doc/2/notebook/notebook.html

From last class:

• we wrote down a budget for the energy content of the global climate system
• we wrote the OLR in terms of an effective emission temperature $T_e$
• The equilibrium emission temperature for Earth is about $T_e \approx 255$ K
• This depends only on energy output from the Sun, and the planetary albedo
• We assumed that global energy content is proportional to surface temperature $T_s$
• We thus have a single equation in two unknown temperatures:

$$C \frac{d T_s}{dt} = (1-\alpha) Q - \sigma T_e^4$$

Parameterizing the dependence of OLR on surface temperature

Later, we will introduce additional physics for column radiative transfer to link $T_s$ and $T_e$. For now, we'll make the simplest assumption we can:

$$ T_e = \beta T_s$$

where $\beta$ is a dimensionless constant. This is a parameterization that we introduce into the model for simplicity. We need a value for our parameter $\beta$, which we will get from observations:

$T_e = 255$ K

$T_s = 288$ K

Thus

$$ \beta = 255 / 288 = 0.885$$

Using this parameterization, we can now write a closed equation for surface temperature:

$$C \frac{d T_s}{dt}=(1-α)Q-σ(\beta T_s )^4$$

Solving the energy balance model

This is a first-order Ordinary Differential Equation (ODE) for $T_s$ as a function of time. It is also our very first climate model.

To solve it (i.e. see how $T_s$ evolves from some specified initial condition) we have two choices:

1. Solve it analytically
2. Solve it numerically

Option 1 will usually not be possible because the equations will typically be too complex and non-linear. This is why computers are our best friends in the world of climate modeling.

HOWEVER it is often useful and instructive to simplify a model down to something that is analytically solvable when possible. Why? Two reasons:

1. Analysis will often yield a deeper understanding of the behavior of the system
2. Gives us a benchmark against which to test the results of our numerical solutions.

Equilibrium solutions

Note that equilibrium surface temperature is

$$\bar{T_s} = \frac{1}{β} \bigg( \frac{(1-α)Q}{σ}\bigg)^{\frac{1}{4}} = 288 K $$

We are going to linearize the equation for small perturbations away from this equilibrium.

Let $T_s = \bar{T_s} + T_s^\prime$ and restrict our solution to $T_s^\prime << \bar{T_s}$. Note this this is not a big restriction! For example, a 10 degree warming or cooling is just $\pm$3.4% of the absolute equilibrium temperature.

Linearizing the governing equation: the feedback parameter

Now use a first-order Taylor series expansion to write

$$OLR \approx \sigma \big(\beta T_s \big)^4 \approx \sigma \big(\beta \bar{T_s} \big)^4 + \Big(4 \sigma \beta^4 \bar{T_s}^3 \Big) T_s^\prime $$

and the budget for the perturbation temperature thus becomes

$$C \frac{d T_s^\prime}{d t} = \lambda T_s^\prime$$

where we define

$$\lambda = -\Big(4 \sigma \beta^4 \bar{T_s}^3 \Big)$$

Here $\lambda$ is actually the feedback parameter for our system. We will say much more about this later in the course. Putting in our observational values, we get

$\lambda = -3.3$ W m$^{-2}$ K$^{-1}$

This is actually our first estimate of what we will later call the Planck feedback, which is found in every climate model. It is the tendency for a warm surface to cool by increased longwave radiation to space. We include the negative sign to indicate that this is a negative feedback, which tends to restore the system towards equilibrium. Again, we will return to this later in the course.

Solve the linear ODE

Now define

$$ \tau = C / (-\lambda)$$

This is a positive constant with dimensions of time (seconds). With these definitions the temperature evolves according to

$$ \frac{d T_s^\prime}{d t} = - \frac{T_s^\prime}{\tau}$$

This is one of the simplest ODEs. Hopefully it looks familiar to most of you. It is the equation for an exponential decay process. We can easily solve for the temperature evolution by integrating from an initial condition $T_s^\prime(0)$:

$$ \int_{T_s^\prime(0)}^{T_s^\prime(t)} \frac{d T_s^\prime}{T_s^\prime} = -\int_0^t \frac{dt}{\tau}$$$$\ln \bigg( \frac{T_s^\prime(t)}{T_s^\prime(0)} \bigg) = -\frac{t}{\tau}$$$$T_s^\prime(t) = T_s^\prime(0) \exp \bigg(-\frac{t}{\tau} \bigg)$$

I hope that the mathematics is straightforward for everyone in this class. If not, go through it carefully and make sure you understand each step.

e-folding time for relaxation of global mean temperature

Our model says that surface temperature will relax toward its equilibrium value over a characteristic time scale $\tau$. This is an e-folding time – the time it takes for the perturbation to decay by a factor 1/e = 0.37

What should this timescale be for the climate system?

To estimate $\tau$ we need a value for the effective heat capacity $C$. A quick and dirty estimate:

$$C = c_w \rho_w H$$

where

$c_w = 4 \times 10^3$ J kg$^{-1}$ $^\circ$C$^{-1}$ is the specific heat of water,

$\rho_w = 10^3$ kg m$^{-3}$ is the density of water, and

$H$ is an effective depth of water that is heated or cooled.

What is the right choice for $H$? That turns out to be an interesting and subtle question. It depends very much on the timescale of the problem

• days?
• years?
• decades?
• millenia?

We will revisit this question later in the course. For now, let’s just assume that $H = 100$ m (a bit deeper than the typical depth of the surface mixed layer in the oceans).

Then $C = 4 \times 10^8$ J m$^{-2}$ K$^{-1}$. And the e-folding time for the surface temperature will be

$$\tau = \frac{4 \times 10^8 J m^{-2} K^{-1}}{3.3 W m^{-2} K^{-1}} = 1.4 \times 10^8 s \approx 4 years $$

This is a rather fast timescale relative to other processes that can affect the planetary energy budget, as we will discuss later.

Some take-away messages:

• Earth (or any planet) has a well-defined equilibrium temperature because of the temperature dependence of the outgoing longwave radiation.
• The system will tend to relax toward its equilibrium temperature on an e-folding timescale that depends on (1) radiative feedback processes, and (2) effective heat capacity.
• In our estimate, this e-folding time is relatively short. In the absence of other processes that can either increase the heat capacity or lower (in absolute value) the feedback parameter, the Earth would never be very far out of energy balance
• We will quantify this statement more as the term progresses.

Plotting the solution in Python

Here I'm going to show some example code for making simple line plots with Python. I strongly encourage you to try this out on your own. Avoid the temptation to copy and paste the code! You won't learn anything that way. Type the code into your own Python session. Experiment with it!

First we'll define a bunch of constants.

This code uses the numpy package to do efficient array operations. Before we use the package, we import it into the current Python session.

To make a plot, we will use the matplotlib library. The plotting commands work a lot like MATLAB. But, like any other package, we need to import it before using it.

Solving the ODE numerically

In this case the equation is sufficiently simple that we have an analytical solution. Most models are too mathematically complex for this, and we need numerical solution methods. Because the governing equations for every model are differential in time (and usually in space as well), we need to use some numerical approximations to the governing equations.

We approximate the time derivative with

$$ \frac{dT}{dt} \approx \frac{∆T}{∆t} = \frac{T_1-T_0}{t_1-t_0} $$

which makes sense as long as the timestep $\delta t$ is sufficiently small.

What is meant by sufficiently small? In practice, small enough that the numerical solution behaves well! We will not spend much time in this course talking about numerical methods, but there is much we could say about this…

The simplest time discretization is called Forward Euler or Explicit Euler. Say we know the state of the system at time $t_0$, i.e. we know the temperature $T_0$. Then rearranging the above,

$$T_1 = T_0 + ∆t (dT/dt)$$

So if we can calculate the tendency of the system (i.e. the time derivative) at time $t_0$, then we have a formula to predict the next state of the system.

For our linearized zero-dimensional energy balance model,

$$\frac{dT_s}{dt} = -\frac{1}{\tau} \big(T_s-\bar{T_s} \big)$$

So we can predict the temperature with

$$ T_1 = T_0 - \frac{\Delta t}{\tau} \big( T_0 - \bar{T_s} \big)$$

Let’s implement this formula as a simple function in Python to calculate the next temperature at each timestep

Now let’s construct the full numerical solution by looping in time:

Now we are going to plot this alongside the analytical solution.

So this works quite well; the two solutions look nearly identical.

Now that we have built some confidence in the numerical method, we can use it to study a slightly more complex system for which we don’t have the analytical solution.

E.g. let’s solve the full non-linear energy balance model:

$$C \frac{dT_s}{dt} = (1-\alpha) Q - \sigma \big(\beta T_s \big)^4 $$

We’ll write a new solver function:

We will now follow the same procedure as above to solve this model and plot the solution along with our two previous solutions of the linear model.

And we see that the models essentially do the same thing.

Now try some different initial conditions

The system relaxes back to 288 K regardless of its initial condition.